%----------------------------------------------------
% Filter for MMG signals
% Edited : 14/01/2013
%----------------------------------------------------
%% Import Initial Data

%D1 = dlmread('D3.txt',',',1,0); %Raw data
Rdata = D3;
l= length(Rdata);
x= Rdata(:,1);
MMG1= Rdata(:,2);
y1 = MMG1;
MMG2= Rdata(:,3);
y2 = MMG2;
data1= y1 - mean (y1); % Normalised data
data2=y2 - mean(y2);

subplot(4,1,1)
plot (x,data1)
xlabel('Time(s)')
ylabel('Amplitude (V)')


%% Filter Details

fs = 1/(mean([x(3)-x(2) x(4012)-x(4011)])/1000000); %Sampling frequency average of two time steps
samplePeriod = 1/fs;
Fc = 2; %cut off frequency
Fc2 = 100;

%% First Order Butterworth High Pass Filter

% wn1 = (Fc)/(1/samplePeriod); %Normalised frequency. RBW's defination
% wn2 =  Fs;
% wn = [wn1 wn2];

w = (Fc)/(1/samplePeriod); %Normalised frequency. RBW's defination

[b, a] = butter(1, w, 'high');

newy =  filtfilt(b, a, data1);
newyy=  filtfilt(b, a, data2);

w = (Fc2)/(1/samplePeriod); %Normalised frequency. RBW's defination

[b, a] = butter(1, w, 'low');

newy2 = filtfilt(b, a, newy);
newyy2 = filtfilt(b, a, newyy);

subplot(4,1,3)
plot (x,data1)
hold on
plot(x,newyy,'r-');
xlabel('Time(s)')
ylabel('Amplitude (V)')

subplot(4,1,4)
plot (x,data2)
hold on
plot(x,newyy2,'r-');
xlabel('Time(s)')
ylabel('Amplitude (V)')

% End
figure;
hold off